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Abstract.  In  this  paper  the  accuracy  and  efficiency  of  a finite- volume 
multigrid  solver  for  Large  Eddy  Simulation  (LES)  is  investigated.  The  spa- 
tial discretization  method  employed  is  a second-order  accurate  central  dif- 
ferencing scheme.  For  time  discretization  of  the  momentum  equations  the 
implicit  second-order  Crank-Nicolson  method  and  the  explicit  second-order 
Adams-Bashforth  method  are  considered.  The  influences  of  the  two  time 
discretizations,  choice  of  grid  size  and  time-step  size  and  multigrid  per- 
formance on  the  numerical  accuracy  and  computational  efficiency  are  dis- 
cussed. 


1.  Introduction 

Due  to  the  foreseeable  progresses  in  the  performance  of  computer  systems, 
it  can  be  expected,  that  in  the  near  future  the  Large  Eddy  Simulation 
(LES)  will  become  more  important  and  applicable  also  in  industrial  prac- 
tice. However,  efficient  numerical  algorithms  designed  for  modern  parallel 
computer  architectures  with  facilities  to  allow  the  modeling  of  complex  ge- 
ometries are  another  crucial  issue  in  the  application  of  LES.  In  this  way  it 
will  be  possible  to  achieve  a sufficient  numerical  resolution  and  geometrical 
flexibility  also  to  deal  with  complex  practical  problems  within  reasonable 
computing  times. 

Numerical  aspects  such  as  spatial  and  temporal  discretization,  solu- 
tion algorithms  and  resolution  requirements  , and  modeling  aspects  such  as 
subgrid  scale  models  are  effecting  the  accuracy  and  efficieny  of  the  simula- 
tions. Due  to  the  enormous  computational  requirements  for  LES,  there  can 
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be  found  few  systematic  investigations  about  the  above  mentioned  aspects. 
Despite  that  different  groups  have  organized  LES  workshops  for  specified 
test  cases,  it  was  difficult  to  come  to  any  definite  conclusions  about  the 
performance  of  different  applied  numerical  methods  as  well  as  the  varying 
subgrid  scale  models. 

In  this  study  the  accuracy  and  performance  of  a finite- volume  multigrid 
solver  for  LES  is  investigated.  The  main  focus  of  the  present  contribu- 
tion will  relate  to  the  following  aspects  with  respect  to  their  influences  on 
the  numerical  accuracy  and  computational  efficiency  of  the  considered  ap- 
proach: comparison  of  the  exp  licit /implicit  methods,  choice  of  grid  size  and 
time-step  size  and  multigrid  performance. 

For  this  an  already  well  investigated  turbulent  channel  flow  will  be  con- 
sidered for  comparison  to  other  results  from  the  literature. 


2.  Numerical  Procedure 


The  governing  equations  for  an  incompressible  flow  are  given  by 


duj 

dxi 

djpuj)  d(pUjUj) 
dt  dxj 
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where  Ui  are  the  velocity  components  of  the  resolved  scales  with  respect 
to  the  Cartesian  coordinates  Xi,  p is  the  corresponding  pressure,  p is  the 
viscosity,  p is  the  density  and  t is  the  time.  The  subgrid  scale  stresses 


Tij  — /?  (ujtij  Ui'Uj'j  (3) 

are  modeled  by  the  Smagorinsky  model  [1]  and  the  dynamic  Germano  model 
[2]  to  close  the  problem.  In  order  to  stabilize  the  dynamic  model,  the  nega- 
tive values  for  the  parameter  Cs  of  the  underlying  Smagorinsky  model  are 
clipped. 

The  basic  flow  solver  is  the  FASTEST-3D  code  (INVENT  Computing, 
Erlangen)  [3]  with  extensions  for  LES.  The  solver  is  based  on  a fully  con- 
servative finite- volume  method  for  solving  the  incompressible  Navier-Stokes 
equations  on  a non-staggered,  cell-centered  grid  arrangement.  The  spatial 
discretization  method  employed  is  a second-order  accurate  central  differ- 
encing scheme  for  block-structured  non-orthogonal  boundary-fitted  grids. 
For  time  discretization  of  the  momentum  equations  we  consider  two  dif- 
ferent approaches:  the  implicit  second-order  Crank-Nicolson  method  and 
the  explicit  second-order  Adams-Bashforth  method.  Both  methods  define 
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approximations  u \ and  f)vh  to  the  solution  of  the  continuous  problem  at  the 
time  levels  tn  = nAt  ( n = 1,2,...),  where  the  parameter  h is  a measure 
for  the  spatial  resolution  and  At  > 0 is  the  time-step  size. 

Within  the  implicit  method,  for  each  time-step,  assuming  that  the  un- 
knowns at  the  time  level  fn_i  have  already  been  computed,  the  unknowns 
at  the  time  level  tn  have  to  be  determined  as  the  solution  of  a nonlinear 
algebraic  system.  For  this  a nonlinear  full  approximation  multigrid  scheme 
with  a pressure-correction  smoother  is  employed  [4], 

The  smoothing  procedure  is  based  on  a variant  of  the  well-known  SIM- 
PLE algorithm  proposed  by  Patankar  and  Spalding  [5].  The  determination 
of  u'l  and  pft  is  done  in  several  steps.  In  the  first  step,  after  the  calculation 
of  the  turbulent  viscosity,  an  intermediate  approximation  to  ufL  is  obtained 
by  solving  the  discrete  momentum  equations  with  the  pressure  term,  the 
source  term  and  the  matrix  coefficients  formed  with  values  of  the  previ- 
ous iteration.  In  the  second  step,  corrections  Ap^  and  Au™  are  sought  to 
obtain  the  new  pressure  pft  and  the  new  velocity  u\  exactly  satisfying  the 
continuity  equation.  By  considering  a modified  discrete  momentum  equa- 
tion together  with  the  discrete  continuity  equation,  an  equation  for  the 
pressure  correction  A is  derived,  where  a selective  interpolation  tech- 
nique is  used  for  making  the  cell  face  velocities  dependent  on  the  nodal 
pressure,  which  is  necessary  to  avoid  oscillatory  solutions  that  may  occur 
owing  to  the  non-staggered  grid  arrangment[6].  To  improve  the  diagonal 
dominance  in  the  pressure-correction  equation  the  contributions  due  to  grid 
non-orthogonality  are  neglected.  The  smoothing  iteration  step  is  completed 
by  correcting  the  velocity  components  and  the  pressure.  To  ensure  conver- 
gence, for  the  velocity  components  and  the  pressure  an  under-relaxation  in 
the  variant  suggested  by  Patankar  [7]  is  employed.  For  the  solution  of  the 
linear  system  of  equations  the  Strongly  Implicit  Procedure  (SIP)  of  Stone 
[8]  is  used.  The  global  outer  multigrid  procedure  is  implemented  as  a non- 
linear full  approximation  scheme,  in  which  the  above  pressure-correction 
scheme  acts  as  the  smoother.  The  different  grid  levels  are  visited  following 
the  standard  V-cycle  approach,  where  second-order  interpolation  is  em- 
ployed for  the  grid  transfers.  This  procedure  is  repeated  till  convergence  is 
reached. 

Within  the  explicit  method,  first  the  turbulent  viscosity  and  the  veloc- 
ity Tifr  are  calculated  explicitly  from  the  unknowns  at  the  time  level  tn-\. 
Then  the  pressure  correction  equation  is  derived  using  the  modified  dis- 
crete momentum  equation  together  with  the  discrete  continuity  equation. 
The  corresponding  resulting  linear  system  is  solved  by  a linear  multigrid 
method.  The  computation  of  the  time  level  tn  is  completed  by  correcting 
the  velocity  components  and  the  pressure.  In  Figure  1 a schematical  flow 
diagram  of  both  methods  is  given. 


3.  Numerical  Results 

The  following  investigations  concern  the  accuracy  and  the  numerical  effi- 
cieny  of  the  methods  given  in  the  previous  section.  All  computations  were 
carried  out  on  a Compaq  AlphaServer  ES40  667  MHz. 

As  a test  case  the  turbulent  flow  between  two  parallel  plates  seperated 
by  a distance  25  is  investigated.  The  flow  is  driven  by  a uniform  streamwise 
pressure  gradient.  The  Reynolds  number  ReT  = 395  (Re\,  = 6875)  based 
on  half- width  and  friction  (or  bulk)  velocity  is  considered.  Since  the  stream- 
wise  and  spanwise  directions  (i.e.  x and  z)  are  formally  infinite,  periodic 
boundary  conditions  are  used  for  the  simulation  of  a finite  domain.  The 
computational  domain  has  the  dimensions  2it5  x 25  x tt5  which  are  con- 
sidered to  be  large  enough  to  avoid  adverse  effects  of  the  periodic  bound- 
ary conditions.  Initial  conditions  with  random  fluctuations  are  given  to 
assure  turbulence.  Three  different  grids  with  64  x 32  x 32,  96  x 48  x 48  and 
128  x 64  x 64  control  volumes  (CVs)  are  employed  in  x,  y,  z directions.  The 
grids  are  equidistant  in  streamwise  and  spanwise  directions  where  in  the 
wall-normal  direction  a geometric  progression  is  used  for  the  grid  density. 
The  first  point  away  from  the  wall  is  at  y+  « 9.25  for  the  coarse  grid, 
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y+  a 6.17  for  the  medium  grid  and  y+  & 4.62  for  the  fine  grid.  Before 
starting  the  computation  of  mean  values  and  statistics  of  the  flow,  the  si- 
mulations are  carried  out  until  the  numerical  solution  reached  a statistically 
steady  state.  All  data  are  computed  from  the  appropriate  quantities  and 
averaged  in  the  homogeneous  spatial  directions  x and  z and  in  time.  The 
averages  in  time  are  taken  over  200  s in  each  case.  The  results  are  compared 
to  the  DNS  data  given  in  [9]. 

First,  a validation  of  the  applied  subgrid  models  is  carried  out.  In  Figure 
2 (left)  a comparison  of  the  mean  velocity  profile  Um  normalized  by  the  fric- 
tion velocity  uT  in  the  normalized  half  channel  width  is  given.  The  results 


Figure  2.  Comparison  of  mean  velocity  profile  (left)  and  rms  velocity  profile  (right) 
between  LES  and  DNS. 


calculated  with  the  finest  grid  and  the  Germano  model  correspond  well  with 
the  reference  data.  Figure  2 (right)  shows  the  root  mean  square  urms  of  the 
velocity  in  streamwise  direction  normalized  by  the  friction  velocity  uT  in 
normalized  half  channel  width  for  three  grids  with  the  Smagorinsky  model. 
Compared  with  the  DNS  data  the  LES  results  for  the  coarse  grid  show 
slightly  lower  fluctuations.  However,  the  overall  agreement  is  satisfactory 
and  with  grid  refinement  the  profile  approaches  the  reference  data. 

Next,  the  numerical  efficiencies  of  the  explicit  and  implicit  methods  are 
investigated.  Figure  3 shows  the  comparison  of  both  methods  without  the 
multigrid  solver  for  CFL=0.5,  where  the  computational  time  per  time-step 
against  the  number  of  CVs  is  plotted.  The  CFL-number  is  defined  by 


CFL  = max 
cv 


■ 3 

^Ui/Axj 

.i=l 


• At 


where  the  maximum  is  taken  over  all  control  volumes  (CV).  It  is  observed 
that  without  a multigrid  solver  the  implicit  method  is  more  efficient  than 
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Figure  3.  Comparison  of  explicit  and  implicit  methods  without  the  multigrid  solver  for 
CFL=0.5  with  varying  grid  size. 


the  explicit  method,  where  the  differences  in  CPU-time  increase  with  grid 
refinement.  This  is  due  to  the  fact,  that  within  the  explicit  method,  the 
pressure-correction  equation  has  to  be  solved  very  exactly,  since  it  is  solved 
only  once.  Thus  the  CPU-time  for  solving  the  linear  system  without  the 
multigrid  solver  is  very  high  and  the  computational  effort  increases  rapidly 
when  the  grid  is  refined. 

In  Figure  4 the  computational  time  per  time-step  against  the  number  of 
CVs  for  both  methods  with  the  multigrid  solver  for  the  same  CFL-number 
is  plotted.  The  multigrid  approach  yields  an  acceleration  for  both  methods. 
The  acceleration  for  the  explicit  method  is  much  higher  as  for  the  implicit 
method,  such  that  the  explicit  method  becomes  superior  to  the  implicit 
one. 


Figure  4 ■ Comparison  of  explicit  and  implicit  methods  with  the  multigrid  solver  for 
CFL=0.5  with  varying  grid  size. 
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Since  there  exists  no  strict  time-step  limitation,  of  course,  with  the 
implicit  method  computations  with  higher  CFL  numbers  are  possible.  In 
Figure  5 the  computing  times  against  the  number  of  CVs  are  plotted  for 
the  implicit  method  using  CFL=6  and  the  explicit  method  using  CFL=0.5. 
It  can  be  seen  that  even  in  this  comparison  the  explicit  method  still  turns 
out  to  be  slightly  more  efficient.  In  order  to  see  how  the  time-step  size  is 
effecting  the  results  of  the  computations,  in  Figure  6 the  root  mean  square 
values  urms  of  the  velocity  in  streamwise  direction  normalized  by  the  fric- 
tion velocity  uT  in  the  normalized  half  channel  width  for  the  coarse  grid 
with  the  Smagorinsky  model  is  shown  for  different  CFL-numbers.  The  com- 
putations with  a larger  time-step  size  (CFL=6)  overpredict  the  streamwise 
fluctuations  when  compared  with  the  smaller  time-step  size  (CFL=0.5). 
Further  test  computations  have  shown  that  using  even  larger  time-step 
sizes  (e.g.  CFL=10)  yield  rather  poor  results,  such  that  for  the  considered 
test  case  CFL=6  can  be  viewed  as  the  maximum  value  to  achieve  physically 
reasonable  results. 


Figure  5.  Comparison  of  explicit  and  implicit  methods  with  the  multigrid  solver  for 
different  CFL-numbers  and  varying  grid  size. 


4.  Conclusions 

The  numerical  accuracy  and  efficieny  of  explicit  and  implicit  methods  with 
and  without  multigrid  for  Large  Eddy  Simulation  has  been  investigated.  For 
the  comparisons  the  simple  well  known  turbulent  channel  flow  is  consid- 
ered. Without  using  the  multigrid  solver,  the  computational  requirement  for 
the  explicit  method  increases  very  rapidly  with  increasing  number  of  grid 
points,  such  that  in  this  case  the  implicit  method  is  superior.  However, 
when  employing  the  multigrid  solver,  a high  acceleration  for  the  explicit 
method  can  be  realized,  which  makes  the  method  more  efficient  than  the 
implicit  one.  It  has  been  shown,  that  with  the  implicit  method  it  is  possible 
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Figure  6.  Comparison  of  rms  velocity  profile  for  different  CFL-numbers. 


to  achieve  also  for  CFL>1  physically  reasonable  results,  but  comparing  the 
overal  computational  performance  the  explicit  method  still  turned  out  to 
be  superior. 

Of  course,  it  would  be  of  interest,  if  the  findings  for  the  considered 
geometrically  very  simple  test  case  will  also  be  confirmed  for  problems  in 
more  complex  geometries.  The  situation  may  change  in  such  cases,  when  it 
is  necessary  to  use  more  irregular  grids  with  higher  variations  of  control- 
volume  sizes.  Here,  the  advantages  of  the  implicit  method,  where  there  is 
no  need  to  adjust  the  time-step  size  to  the  smallest  control  volumes,  can 
become  more  dominant.  Corresponding  investigations  will  be  a topic  of 
forthcoming  work. 
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